Critical Behavior of Ferromagnetic Ising Model on Triangular Lattice 
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We apply a new updating algorithm scheme to investigate the critical behavior of the two- 
dimensional ferromagnetic Ising model on a triangular lattice with nearest neighbour interac- 
tions. The transition is examined by generating accurate data for large lattices with L — 
8, 10, 12, 15, 20, 25, 30, 40, 50. The spin updating algorithm we employ has the advantages of both 
metropolis and single-update methods. Our study indicates that the transition to be continuous 
at Tc — 3.6403(2). A convincing finite-size scaling analysis of the model yield v — 0.9995(21), 
P/u = 0.12400(18), = 1.75223(22), -y' /u = 1.7555(22), a/u = 0.00077(420) (scaling) and 

a/v = 0.0010(42) (hyperscaUng) respectively. Estimates of present scheme yield accurate estimates 
for all critical exponents than those obtained with Monte Carlo methods and show an excellent 
agreement with their well-established predicted values. 

PACS numbers: 75.10.Nr, 64.60.Fr, 75.40.Cx 



I. INTRODUCTION 

Ising Model on triangular lattice, as an archetypical 
example of a frustrated system, was initially studied by 
Wannier [l[ and NewU The triangular Ising model 
has attracted much attention and because of its remark- 
able properties it has a long history of investigation. No 
exact solution is available in two dimensions in an ar- 
bitrary magnetic field. Hence, simulations of the Ising 
model are essential. Monte Carlo simulation methods 
have been widely using techniques to update the spins of 
the system and to study the Ising model on triangular 
lattice to obtain numerical solutions. 

A number of Monte Carlo methods based on Metropo- 
lis algorithms 0] have been applied to the model in the 
past with somewhat mixed results. A classical Monte 
Carlo using Metropolis algorithm runs into difficulties on 
large lattice sizes due to critical slowing down and rapid 
increase in correlation time. Because of these reasons 
it is difficult to obtain meaningful results on larger lat- 
tice sizes. A cluster technique, which uses multi-cluster 
algorithm, was pioneered by Swendsen and Wang [J]. 
Their method demonstrated its validity and efficiency 
for the Potts model successfully. Similar ideas have been 
pursued by Wolff Q who proposed single-cluster algo- 
rithm, a nonlocal updating technique based on multiple- 
cluster algorithm but more efficient and easily applicable 
to achieve the meaningful results. Both of these cluster 
algorithms are very effective in reducing critical slowing 
down. But for the same system, one sweep of single- 
cluster or multiple-clusters take much more time than 
that of Metropolis algorithm, hence less effective for large 
lattice sizes. Very little has been done since then using 
these approaches on this model. Since large systems have 
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the advantage of reducing finite-size effects, for this rea- 
son, we are forced to look yet again for an alternative 
approach. 

In this study, we attempt to use a mixture of extend 
a single-cluster (multiple-cluster) and Metropolis algo- 
rithms to the Ising model in two-dimensions on trian- 
gular lattice. Applications of combined algorithm tech- 
niques have been extremely successful in lattice QCD 
[10, i, B [13, [m, [H, [11 and have given rise to great 
optimism about the possibility of obtaining results rele- 
vant to continuum physics from Monte Carlo simulations 
of lattice version of the corresponding theory. As men- 
tioned above, our aim is to use standard Monte Carlo 
techniques based on combined effect of single-cluster 
(multiple-cluster) and Metropolis algorithm to update 
the spins of the system and see whether useful results 
can be obtained. The values of the critical exponents 
are well established for this model in 2-dimensions and 
in order for our method to be considered successful, it 
must reproduce these well-established values. For this 
purpose, the critical exponents are computed by using 
the combined algorithm. 

The rest of the paper is organised as follows: In Sec. II 
we discuss the Ising model in 2-dimensions in its lattice 
formulation. Here we describe the details of simulations 
and the methods used to extract the observables. We 
present and discuss our results in Sec. HI. Our conclu- 
sions are given in Sec. IV 

II. MODEL AND SIMULATION DETAILS 

The Hamiltonian of the model is given by 

h = -jY, S,S,, (1) 

<l,j> 

where J is positive and denotes the strength of the ferro- 
magnetic interaction, Si = ±1 is the Ising spin variable. 



2 



and < i,j > restricts the summation to distinct pairs of 
nearest neighbors. The term on the right-hand side of Eq. 
[T] shows that the overall energy is lowered when neigh- 
bouring spins are aligned. This effect is mostly due to 
the Pauli exclusion principle. No such restriction applies 
if the spins are anti-parallel. The geometrical structure 
of the model is shown in Fig. [T] Periodic boundary con- 
ditions were employed during the simulations to reduce 
finite-size effects. 



FIG. 1: Triangular Ising net. 



Configuration ensembles were generated using both 
mixture of Metropolis and single-cluster (multiple- 
cluster) methods on two dimensional triangular lattice 
of size L = 8 — 50. We outline the procedure of single- 
cluster updating algorithm used to study the Ising model 
on triangular lattice below: 

(a) Choose a random lattice site x as the first point of 
cluster c to built. 

(b) Visit all sites directly connecting to x, and add 
them to the cluster c with probability P{Sx,Sy) = 
1 — exp{min[0,2KSxSy]}, where y is the nearest 
neighbor of site x and K is the inverse tempera- 
ture. 

(c) Continue iterativcly in the same way for the newly 
adjoined sites until the process stops. 

(d) Flip the all sites of the cluster c. 

(e) Repeat(a), (b), (c) and (d) sufficient times. 

We define a compound sweep as one single-cluster up- 
dating sweep following with five Metropolis updating 
sweeps. In our simulations, five compound sweeps were 
performed between measurements. We observed that 
Metropolis - single-cluster (multiple-cluster) and com- 
pound sweep updating techniques proved equally good 
for large ensembles. However, in case of the compound 
sweep technique the computational cost turned out to be 
far less. For this reason we adapted this technique in our 
simulations. 

We generated 1 x 10^ ensembles at transition point 
for each lattice size to performed the finite-size analysis. 
High statistics could help us to obtain a more accurate 
result. Since the ensembles were generated by Monte 



Carlo method, they were not completely uncorrelated. 
The Jackknife procedure was employed because it took 
the autocorrelation between the ensembles into account 
and provided an improved estimate of the mean values 
and errors. 



III. RESULTS AND DISCUSSION 

From the fluctuation dissipation theorem, the heat ca- 
pacity (C) is the derivative of the energy with respect to 
temperature and has the form 



C 



< E'^ > - < E>' 
L2T2 



(2) 



For improved comparability among different system sizes, 
it is better to compute the heat capacity per spin which 
is displayed in in Fig. [21 Theoretical derivations suggest 
that the heat capacity should behave like log \ T — \ 
near the critical temperature. We notice a progressive 
steepening of the peak as the lattice size increases, which 
illustrates a more apparent phase transition. Note that 
due to finite size effects, the peak is flattened and moved 
to the right. Our interest is to look for such a transition 
point and investigate the critical phenomena associated 
with the transition using finite-size scaling method [3, 

MM 




temperature 



FIG. 2: Plot showing the differing results of the heat capacity 
with respect to temperature for varying lattice size, L x L. 



A. Estimation of the transition temperatures 

The fourth-order magnetic cumulant Ul , which is used 
to estimate the transition temperature, is defined by the 
following expression [l6| : 



f/(r,L) = i 



< m^{T,L) > 
3 < m^{T,L) >2 



(3) 
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where < m'^ (T, L) > is the thermodynamic average value 
of the kth power of the magnetic order parameter per spin 
for the lattice of size L with temperature T. The varia- 
tion of U (T, L), for a given size L, with the temperature 
is illustrated in Fig. [3l To determinate the temperature 
Tc{L, L'), we make use of the condition [H, [l3| 



U{T,L') 



U{T, L) 



1, 



(4) 



T=T^{L,L' 



where L' and L are two different lattice sizes. Thus we 
can determine Tc{L,L') by locating the interaction of 
these curves - the so called "cumulant crossing" . 
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FIG. 3: The fourth-order magnetic parameter cumulant plot- 
ted as a function of inverse temperature for several lattice 
sizes. 

Fig. m shows the results for Tc{L,L') plotted as a 
function of the ratio ln~^(L'/L). Keeping L fixed, lin- 
ear extrapolation is performed to obtain Tc{L, oo) which 
corresponds to the critical temperature of lattice size L. 
Results of the extrapolations for 8 < L < 15 agree quite 
well within the errors. The transition temperature for 
the infinite lattice is thus estimated as Tc = 3.6403(2), 
which is much close to the exact value 4/ In 3. 



B. Order of the transition 

A better way to identify the order of the transition 
is the internal-energy cumulant [H, [H, HO, HH [H [H 
defined by 



V{T,L)^\- 



< E^{T,L) > 
3 < E^{T,L) >2 



(5) 



where < E^{T, L) > is the thermodynamic average value 
of A:th power of internal energy for the lattice of size L 
with temperature T. V(T, L) is a useful quantity since 
its behavior at a continue phase transition is quite dif- 
ferent from that at a first order transition. This variable 




FIG. 4: Estimates for Tc plotted vs inverse logarithm of the 
scale factor b = L' /L for serval lattice sizes. 



has a minimum (near the critical point) , and in the limit 
L — > oo achieves the value ^* = f for a continuous 
transition, whereas V* < ^ is expected in the case of a 
first-order transition. Plots of V(T, L) versus tempera- 
ture for several lattice sizes are shown in Fig. [5l The 
graphs indicate that the transition is not first-order since 
the energy culumant V{T, L) does not peak near the crit- 
ical temperature. This has also been observed in various 
other studies [13, [13, HI . 

Fig. [6] shows the order parameter cumulant U{T,L) 
as a function of temperature for serval lattice sizes. In 
contrast with obtaining a negative minimum value in first 
order transition, U{T,L) drops from 2/3 for T < Tc to 
for T > This is in agreement with the continuous 
transition UM. 
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FIG. 5; Plots of energy cumulant with respect to temperature 
for serval lattice sizes. 

Since the magnetic susceptibility x scales as L'^ for 
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FIG. 6: plots of order parameter cumulant with respect to 
temperature for serval lattice sizes. 



a discontinuous transition and as L''^'^ for a continuous 
one [25j . we can also determine the transition order by 
x's scaling behaviors. As discussed below, we find that 
X scales as Li-75223 ^j. t <Tc and as ^j. y > j^^^ 

but not as L^. Thus wc can conclude that the transition 
is continuous. 



C. Estimation of the critical exponents 



plot of the scaling relation corresponding to this quan- 
tity gives the correlation length exponent ly ~ 0.9995(21), 
which is slightly different from the theoretical value. 
Since the critical properties of ln{dU/dK) have a large 
lattice size dependence [l^, we could also extract the 
value of v from other observables, which are defined as 

following [H,!!!,!!!]: 



Vi 



In 



din < m > 



dK 



In 



<E> 



< mE > 



< m > 





ln( 


E{< Tn> -r 




< m > 


V2 = 


In 


din < 771^ > 
dK 


Vs = 


2[m] — [to^]. 


Vi = 


3[m2] - 2[m% 


where 







(8) 



^^^ Ei<m^>-m^) ^^ (9) 



< TO^ > 



(10) 

(11) 



In 



d < to" > 
dK 



ln(< S >< to" > - < to"£: >) 
In < E{< to" > -to") > . 



(12) 



The curves show a smooth behavior and linear fitting of 
the data is adopted. We found the slopes of the straight 
lines of Vi, with i ~ 1,2,3,4, are very close to that of 
\n{dU/dK), as illustrated in Fig. [71 



To get an estimation for the critical exponents, finite- 
size scaling relations at critical point arc used. For this 
purpose, we extract the critical exponents, which are re- 
lated to specific heat, the order parameter, and the sus- 
ceptibility, from our simulation data at Tc(oo) by using 
finite-size analysis [ilEIli. 

First we extract i/ from dU /dK which obey the follow- 
ing relation 



dK ^ ' 



(6) 



Finite-difference derivative method can be used to deter- 
minate dU / dK, but it is rather a poor choice since the 
results are highly sensitively to the interval of the linear 
approximation. To obtain a precise result of the criti- 
cal exponent v, one can write the derivative of U with 
respect to K as 




FIG. 7: The critical exponent v is determined from the slope 
of \TL{dUL/dK) and Vi vs InL at T = 3.6403. Linear fittings 
result in = 0.9995(21). 



dU 
dK 



1 



3 < TO^ >2 

+ < ni^E >] 

m^i< E > +E) 

< ^ > 

3 < >2 



< TO >< E > -2 



< to"* >< m^E > 
< mP > 



2ni^E < > , , 
3 < 



Figure [5] displays the double-logarithmic plot of the 
magnetization at transition point as a function of L. Ac- 
cording to the standard theory of finite-size scaling, at 
Tc the magnetization per spin should obey the relation 



(13) 



The plot of ln(d[/ jdK) versus In L is shown in Figure [T] for sufficiently large lattice size. The data falls nicely on a 
The slope of the straight line obtained from the log-log linear curve and the slope of the curve gives an estimate 
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for the ratio of critical exponents j3/v = 0.12400(18). 
Using our estimate oi v ~ 0.9995(21), we get the critical 
exponent of magnetization j3 ~ 0.12394(32). 




FIG. 8: Plot of In m as a function of InL at T = 3.6403. 

On finite lattices, the magnetic susceptibility per spin 
is given bv[T7l.[27t 



for T>Tc and 

X' = {LVT){< 
for T < Tc, and satisfy 



> 



(14) 



(15) 



(16) 



As another estimation, we compute the magnetic sus- 
ceptibility per spin by examining log-log plot shown in 
Fig. [9] of X a-nd x' versus L at the transition point. 
The data displays a smooth scaling behaviour. The 
slope of the linear fit to the data gives estimates of 
= 1.75223(22) and y/iy = 1.7555(22), respectively. 
With our estimated values = 0.9995(21), we obtain 
7 = 1.7514(37) and 7' = 1.7546(43), respectively. These 
results agree rather well with the theoretical prediction. 

With three critical exponents v, 7, /3 determined, the 
fourth exponent a is estimated using scaling law 



2/3- 



-7 = 2. 



(17) 



This yields a/j^ = 0.00077(420) and a = 0.00077(420) 
for j/v = 1.75222(22), a/iy = -0.0025(47) and a = 
-0.0025(47) for j'/iy = 1.7555(22) respectively. We can 
also estimate a from hyperscaling relation 



diy 



(18) 



where d = 2 is dimension of the system. It gives a result 
a = 0.0010(42) and a/v = 0.0010(42), which is consis- 
tent with that given by scaling law. 

Estimates of critical exponents from finite-size scaling 
at critical temperature Tc are summcrized in Table ID 




FIG. 9: Plot of Inx and Inx' vs InL at T = 3.6403. 



TABLE I: Estimates of critical exponents 



exponent 


value 


exact value [28. 291 


V 


0.9995(21) 


1 


P/v 


0.12400(18) 


1/8 


l/v 


1.75223(22)(r > Tc) 


7/4 


ilv 


1.7555(22)(r < Tc) 


7/4 


ctjv 


0.00077(420) (by scaling) 





ajv 


0.0010(42) (by hyperscaling) 






To ensure the validity of our estimate of a, we com- 
pare it with those obtained using the specific heat per 
spin from the fluctuations of the total energy. For a con- 
tinuous transition it behaves as 



C 



bL"/". 



(19) 



Using ajv ~ 0.00077(420) calculated above, the heat 
capacity at critical temperature is plotted as a function of 
^0.00077 --^^ Y\g. [ini It was found that the evaluations of 
(fTT)) and (fig]) yield results very consistent within statisti- 
cal errors. We find the values are in good agreement with 
the theoretical ones. An over all error of about 0.1 — 0.8% 
is estimated for the values of the critical exponents 



IV. CONCLUSIONS 

We have obtained critical exponents for the Ising 
model with nearest-neighbour interactions on a trian- 
gular lattice by using a mixture of single-cluster and 
Metropolis algorithm in our simulations. The data are 
analysed according to the finite-size scaling theory. The 
idea of using extensions of mixed algorithm to estimate 
the critical exponents seems to supply a quite accurate 
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FIG. 10: Plot of heat capacity C with respect to L"'''', where 
a/v = 0.00077 and C is evaluated at Tc = 3.6403. 



route for their estimation. In conclusion, it can be stated 
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